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The capability of density-functional theory to deal with the ground-state of strongly correlated 
low-dimensional systems, such as semiconductor quantum dots, depends on the accuracy of function- 
als developed for the exchange and correlation energies. Here we extend a successful approximation 
for the correlation energy of the three dimensional inhomogeneous electron gas, originally intro- 
duced by Becke [J. Chem. Phys. 88, 1053 (1988)], to the two-dimensional case. The approach 
aims to non-empirical modeling of the correlation-hole functions satisfying a set of exact properties. 
Furthermore, the electron current and spin are explicitly taken into account. As a result, good 
performance is obtained in comparison with numerically exact data for quantum dots with varying 
external magnetic field, and for the homogeneous two-dimensional electron gas, respectively. 

PACS numbers: 73.21. La, 71. 15. Mb 



I. INTRODUCTION 

Density-functional theory 1 (DFT) maps the compli- 
cated many-particle problem onto a simple one of non- 
interacting electrons moving in an effective, local poten- 
tial, the Kohn-Sham (KS) potential. The latter is con- 
structed in such a way that the ground-state density of 
the non-interacting particles reproduces the ground-state 
density of the interacting system. Practical success of 
the approach depends on finding good approximations 
for the exchange-correlation (xc) energy which, through 
functional derivation with respect to the particle den- 
sity, defines the xc part of the KS potential. Most of the 
approximations developed so far have focused on three- 
dimensional (3D) systems, i.e., atoms, molecules, and 
solids. Many advances have been made beyond the com- 
monly used local (spin) density approximation [L(S)DA] 
by means of, e.g., generalized gradient approximations, 
orbital functionals, and hybrid functionals. 2 Such efforts 
for two-dimensional (2D) systems have been relatively 
scarce despite the rapidly increasing experimental and 
theoretical interest in quasi-2D structures such as semi- 
conductor layers and surfaces, quantum-Hall systems, 
graphene, and various types of quantum dots 3 (QDs). 

When using DFT, QDs are most commonly treated 
using the 2D-LSDA exchange 4 combined with the 2D- 
LSDA correlation parametrized first by Tanatar and 
Ceperley 5 and later, with more satisfactory spin depen- 
dence, by Attaccalite et al. 6 In many cases, the LSDA 
(prefix "2D" omitted below) performs relatively well 
compared, e.g., with quantum Monte Carlo calculations. 7 
Nevertheless, there is a lack of 2D functionals to deal with 
diverse few-electron QD systems, especially in the strong- 
correlation regime. Only very recently, a local correlation 
functional was developed in 2D within the Colle-Salvetti 
approach, which was found to outperform the LSDA. 8 
In its current form, however, this local functional ap- 



plies only to closed-shell systems with zero spin and zero 
current. In addition to the correlation-energy functional, 
new exchange- energy functionals have been developed for 
finite 2D systems in our foregoing works. 9-11 

In this paper we develop a correlation-energy func- 
tional in 2D. In the derivation, along the lines of the 
work of Becke 12 for the 3D case, we introduce a model 
for spin-dependent correlation-hole functions satisfying a 
set of exact properties in 2D. As a result, we find a spin- 
and current-dependent approximation for the correlation 
energy. In comparison with numerically exact results, the 
obtained accuracy is found to be superior to the LSDA. 
However, we also find, and elucidate, that further model- 
ing of the dependency on the average electron density of 
those parameters describing the size of the correlation- 
hole functions in terms of the size of the exchange-hole 
(x-holc) functions would be required. The applications 
of the functional to a set of few-electron QDs with var- 
ious relative amounts of correlation, ground-state spins, 
electron currents, and external magnetic fields, confirm 
the overall usefulness of the approach. 

II. THEORY 

Within the Kohn-Sham (KS) method of spin-DFT, 13 
the ground-state energy and spin densities pj (r) and 
Pl (r) of a system of N = N-\ + N± interacting electrons 
are determined. The total energy of the interacting sys- 
tem is written as a functional of the spin densities 14 

E v [pt,Pl\ = T s[PhPl\+ J drv(r)p(r) 

+ E H [p}+E xc [p hPl } (1) 

where T s [p-\,pi\ is the kinetic energy functional of non- 
interacting electrons with spin densities p\, p^. v is 
(at vanishing external magnetic field) the external (lo- 
cal) scalar potential acting upon the interacting system, 
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Eh[p\ is the classical electrostatic or Hartree energy of 
the total charge density p = p-f + Pj., and E xc [pf,p±] is 
the xc energy functional. E xc may be further decomposed 
into the exchange energy, E x , and correlation energy E c . 
We have already considered E x in Ref. 9, and thus we 
here focus on E c . Our starting point is the formal ex- 
pression for E c in terms of the correlation-hole (c-holc) 
function, 
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h%° ( ri ,r 2 )= / dA^(n,r 2 ) 
Jo 

= / rfA/ir'(ri,r 2 )-/iS(ri,r 2 )^ (3) 
Jo 

is the c-hole function. Here {aa 1 } = {tt, ii> ti, it} 
and the parameter A G [0, 1] is the electronic coupling 
strength. 1 In the above expression, the c-hole corresponds 
to the full hole function h™ subtracted by the x-hole 
function defined as 

K(r u r 2 ) = t—r— . (4) 

Ax( r 

Note that here h% is defined to have a negative sign in 
contrast to the standard definition [see, e.g., Eq. (4) in 
Ref. 9]. The A-dependent hole function is given by 



. P%% (ri,r 2 ) 
K" (ri,r 2 ) = — -r p CT '(r 2 , 
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where P^J is the two-body reduced density matrix 



(2BRDM) 
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P?x (n,r 2 ) 



N(N -1) dS dA 



dN 



2, JV)* A (1, 2, JV). (6) 



Here, ^(l, 2, N) stands for the ground-state many- 
body wavefunction which is the exact solution of the elec- 
tronic system with a Coulomb coupling strength A, J dN 
denotes the spatial integration and spin summation over 
the iV-th spatial spin coordinates (rjv, ctjv), and we have 
identified a\ = <r, <r 2 = a'. Note that the spin densities 
in Eq. (5) are the same as in the actual, fully interacting 
system. 

As it is well known, in determining the correlation en- 
ergy it is sufficient to know the angular average of the 
c-hole function. In 2D is natural to consider the cylin- 
drical average, given by 



(r. 



i r 27r 
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dep < CT A (r,r + s), 



(7) 



where ri = r, r 2 = r + s, and <j> is the angle between r 
and s. The modeling of the c-hole functions is based on 
satisfying the following conditions: 



sum rule 



(r,*)=0; 



(8) 



• correct short-range behavior for s — > 0; 

• a proper decay in the limit s — > oo; 

• characteristic size assumed to be proportional to 
the characteristic size of the corresponding angular 
averaged x-hole. 

The short-range behavior of the c-hole can be worked 
out by considering the electronic cusp conditions for the 
2D electronic wave function. 15 The angular average of 
the 2BRDM in Eq. (6) is then given by 



1 + I XS 



and 



0) «A CTff (r)(l + 2As), 



(9) 



(10) 



for the same-spin and opposite-spin elements, respec- 
tively. Note that the symbol a always denotes spin oppo- 
site to <7, whereas a', when used, can be equal to cither 
a or a. Here the coefficients A aa i(v) have a spatial de- 
pendence. Although the cusp condition as developed in 
Ref. 15 pertains to the center-of-mass and relative coor- 
dinate, it can be shown, as noted by Becke, 12 that, at 
the order considered above, the same results apply also 
for the coordinate system of r and s. 

The short-range behavior of the 2D x-hole is known. 9 ' 16 
Thus, it is possible to consider the following model for 
the same-spin and opposite-spin c-hole functions, respec- 
tively: 



x F( 7(TCT (r)s), 



B aa {v) - D a (r) + -\B aa {Y)s 



(11) 



and 



K%{v,s) = [B as (r)~ p s {v) + 2\B aS (v) S ] 

x F{ las {v)s). (12) 

Here B aa := A aa j p a and B aS := A aS j p a arc coefficients 
to be determined, and 



_if i (Vp.) 2 i 2 P ,o 

Ua -2\ T ° 4 p^ Pa 



(13) 



where r a 



Ef=i|V^fc, CT | 2 is the (double 



of the) kinetic-energy density, and 



is 



Jp,cr — 

the spin- 
dependent paramagnetic current density. In Eqs. (11) 
and (12), the functions F(~f aa > (r) s) are introduced to 
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ensure the decay of the c-holes in the limit s — > oo. We 
choose them to have the form F(x) = cxp(— a; 2 ), which 
seems appropriate in the case of finite 2D systems. The 
parameters 7 CTCT '(r) are determined by the zero-integral 
constraint of Eq. (8) (see below). 

Next, we introduce the characteristic sizes of the c- 
holes, z aa and z aS , for which the corresponding c-hole 
function vanishes, i.e., 



and 



h a ° x {T,Z aa )=Q\ 



(14) 



Implying these conditions, the coefficients B aa and B a 
in Eqs. (11) and (12) can be written as 



B„ 



1 + §Az 



B n 



1 + 2\z a 



(15) 



It should be noted, however, that the size parameters z aa 
and z aS are functions of r. As for the 3D functional, we 
assume the size of the c-hole to be proportional to the 
size of the x-hole. Thus, we set 

= 2c aa \U° x (v)\-\ (16) 



and 



Zas {r) :=c a - a [ICTOr 1 + |E£(r) 



(17) 



where U£ is the x-hole potential 9 ' 17 for spin a, and c aa 
and c aS are constants to be determined (see below). As 
suggested by Bcckc, 12 the proportionality is plausible due 
to the fact that each electron is surrounded by its Fermi 
hole, and hence the electrostatic interaction between two 
electrons can be expected to be screened beyond some 
characteristic length proportional to the average size of 
the x-hole. The argument is not spin-related, beyond the 
fact that the characteristic length for the x-hole could 
be different for spin up than for spin down. Of course, 
this is only a simplification, but the physical picture is 
appealing, and it has led to good results in 3D atomic 
systems. 

The A-dependent c-hole functions can now be written 

as 



K%V,8) = -\D a {v) 



r (r) 



and 



(r, s) = 2\p- a (v) 



1 + §Az CTCT (r) 



s - z aS (v) 



s 2 F{ laa {v) a), 
(18) 

F( laS (r) S ). (19) 



l + 2\ Z<jS (r) 
Integrating over A yields 

Ax(r) [s-^ CT (r)]s 2 F( 7CTa (r)s) 



K°{r,s) = 



2 Z(j (j (r) 



2^(r) 



3 In [ -z ffcr 



(r) + l 



(20) 



hf(r,s) 



p s (v) [s- ZijS {y)]F (y(r)s) 
2^(r) 



x [2^(r)-ln(2^(r) + l)] . (21) 
Finally, we enforce the sum rules in Eq. (8) giving 

30? 



7<ra(r) 



4^(r) ! 



and 



7<ra(r) = 



2z ffff (r)' 



(22) 



(23) 



This concludes the derivation of the c-hole functions h?" 
and h"" , apart from the determination of constants c aa , 
and c CT ct! respectively. 

From the c-hole functions we can calculate the c-hole 
potentials as 



/•OO 

Uf{v) = 2n dshf{v,s). 
Jo 



(24) 



For the same- and opposite-spin cases of our approxima- 
tion, we find respectively 



Ifi 

Ur(r) = — (8-37r)^(r)z 2 CT (r) 

Ol7T 



2z CTCT (r)-31n( -z aa (r) + 1 



(25) 



and 



Uf(r) = (2-n) Ps (r) 

x [2^(r)-ln(2^(r) + l)]. 

The correlation energies are given by 

E f = \ J dv Pa {v)Uf{v). 



Thus 



E c [p hPl ]=Er +Et L +2E ( 



(26) 



(27) 



(28) 



where we have used the condition = E^ . 

Alternatively, we can compute the correlation energy 
directly from the c-hole functions. From Eqs. (2) and 
(7) we get 

E c [p hPl ]=nJ2 j dv J ds Pa {v)hf{v,s). (29) 

(7<t' 

We remind that D a (r) introduced in Eq. (13) vanishes 
for all the single-particle (N = 1) systems. 12 There- 
fore, the c-hole and thus E c vanish as well. In other 
words, our approximation for the correlation energy is 
self- interaction free for N = 1. 26 
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TABLE I: Total energies from the full configuration- 
interaction calculations (Ref. 21) for totally spin-polarized 
(S z = N/2) quantum dots, exact-exchange total energies, the 
reference correlation energies (El = E tot — Ef^ yi ), c aa yield- 
ing E™ odcl = E T c ci , E™ odcl obtained with a fixed average value 
c aa = 1.32, and the LSDA correlation energy. 



N 


UJ 


Ijircf 
^tot 






Cera 


rnmodcl 

E c 


£,LSDA 


3 


1/4 


2.081 


2.103 


-0.0226 


1.27 


-0.0245 


-0.0538 


3 


1/16 


0.6908 


0.7075 


-0.0167 


1.41 


-0.0144 


-0.0382 


6 


1/4 


7.233 


7.296 


-0.0640 


1.24 


-0.0750 


-0.1125 


6 


1/16 


2.553 


2.599 


-0.0458 


1.36 


-0.0441 


-0.0795 



III. NUMERICAL RESULTS 



TABLE II: Similar to Table I but for unpolarized (S z = 0) 
quantum dots. The correlation energies from our functional, 
_E™ odel , have been calculated using the fixed average values 





= 1.32 and c aa = 


0.75. 










N 


UJ 


zriref 
-C-tot 


r,EXX 


Ef l 




j-Tixiodcl 

E c 


t-iLSDA 
E c 


2 


1 


3* 


3.162 


-0.162 


0.72 


-0.171 


-0.199 


2 


1/4 


0.9324 f 


1.046 


-0.114 


0.82 


-0.102 


-0.139 


2 


1/16 


0.3031+ 


0.373 


-0.070 


0.96 


-0.053 


-0.085 


6 


1/4 


6.995+ 


7.391 


-0.396 


0.73 


-0.406 


-0.457 


12 


1/1.89 2 


25.636* 


26.553 


-0.917 


0.71 


-0.983 


-1.000 



* Analytic solution by Taut from Ref. 18. 

f CI data from Ref. 21. 

t Diffusion QMC data from Ref. 23. 



The first task in the numerical applications is to com- 
plete the correlation functional by finding approxima- 
tions for constants c aa and c a s in Eqs. (16) and (17). 
For this purpose, we consider a set of harmonically con- 
fined QDs, where the external confinement is given by 



v(r) 



,2„2 



/2. Reference results for the correlation en- 



ergies can be obtained from 

E„ — E in , 



E. 



EXX 



(30) 



where El° l t is the exact total energy obtained, e.g., from, 
an analytic, accurate configuration-interaction (CI) or 
quantum Monte Carlo (QMC) calculation, and EXX 
refers to the exact exchange. Here we have calculated 
the EXX energies in the Krieger-Li-Iafrate 19 (KLI) ap- 
proach 19 in the octopus DFT code. 20 The self-consistent 
EXX result - the x-hole potential, (spin) density, kinetic- 
energy density, and current density - is used as input for 
our correlation functional. 

Table I shows the results for the same-spin case. Now 
the QDs are completely spin-polarized with S z = N/2. 
For each QD we show that value of c aa which yields the 
reference correlation energy, i.e., E™ odcl = E r c ci . In these 
examples we find c aa — 1.24 . . . 1.41. Thus, the variation 
of c aa is rather small in view of the fact that the density 
parameter, defined in harmonic QDs as r s = A^~ 1//6 aj~ 2 / 3 
(Ref. 22), varies from 1.9 to 5.3. The second last column 
of Table I shows the correlation energy obtained by using 
a fixed average value of c aa = 1.32. This leads to the 
maximum deviation of ~ 23 % from E 1 ^ 1 . In comparison, 
the self-consistent LSDA correlation energy (last column) 
deviates from the reference result by up to 130 %. 

Table II shows the results for a set of unpolarized 
(S z = 0) QDs in the range 1.5 < r s < 5.7. Note that 
for N > 2 both same- and opposite-spin components of 
the correlation are present, and we have fixed c aa = 1.32 
according to the conclusions above. Fixing c aS = 0.75 
yields deviations of only < 10% from E r c ai , except for 
the highly correlated case of N = 2 and uj = 1/16, which 
shows a deviation of 25 %. The LSDA is still considerably 
further off the reference result than our functional, but 
it performs relatively better than in the polarized case 
discussed above. In particular, when N = 12 the error in 
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FIG. 1: (color online). Total energies (minus the confine- 
ment energy) for a six-electron quantum dot as a function of 
the magnetic field. Results are shown for the reference quan- 
tum Monte Carlo calculations from Ref. 7 (red solid line), for 
the exact-exchange data (blue dashed line), for our functional 
with Cera = 1-32 and c aa = 0.75 (circles), for our functional 
with Caa = 1.1 and Caa = 0.7 (crosses), and for the local-spin- 
density approximation (green dotted line). The total energy 



is calculated from our functional as EF, 



= E u 



+ E" C 



the LSDA correlation is only about 9 %. This is in line 
with the well-known fact that both the L(S)DA exchange 
and correlation, respectively, become more accurate with 
increasing particle number. 

In Fig. 1 we consider the total energy of a more gen- 
eral case: A six-electron QD as a function of the mag- 
netic field B directed perpendicular to the dot plane. 
Increasing the field leads to non-trivial changes in the 
ground-state quantum numbers (S Z ,L Z ) and hence to 
"kinks" in the ground-state total energy as a function 
of B. As the reference data, we use here the variational 
QMC results (red solid line) given in Refs. 7 and 25 for 
a wide range of B up to total spin polarization. The 
confinement strength is here uj = 0.42168, correspond- 
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ing to a typical confinement of 5 meV when modeling 
QDs in GaAs. 3 Note that the total confinement energy, 
i.e., 617 = 6-y/w 2 + where uj c = B/c, has been sub- 

tracted from the total energies to clarify the comparison. 
We point out that the variational QMC method gives 
an upper bound for the true total energy. On the ba- 
sis of previous comparisons between the variational and 
diffusion QMC, and exact diagonalization, 24 our refer- 
ence data in Fig. 1 can be expected to overestimate the 
exact total energy by at most 0.2 . . . 0.3 meV. The max- 
imum possible errors are smaller in the polarized regime 
(B > 5) T. 

Overall, Fig. 1 shows reasonable agreement of our func- 
tional (circles) with the QMC data through the full range 
of the magnetic field. However, the functional yields sys- 
tematically slightly too low correlation energies, and thus 
too low total energies, even if the possible overestima- 
tion of the total energy given by the variational QMC is 
taken into account. We point out that the functional is 
here applied with the fixed parameters c aa = 1.32 and 
c aS = 0.75 suggested by the results in Tables I and II, 
which correspond to considerably weaker confining po- 
tentials (smaller values of to). Obviously, this difference 
implies a different average electron density, and thus a 
different range of the relative correlation energy. In fact, 
if the parameter values are reduced to c a a = 1-1 and 
c aS = 0.7, excellent agreement with QMC is found (see 
the crosses in Fig. 1). Hence, it seems that particularly 
high precision of our functional would require model- 
ing of c aa and c a9 as a function of the particle density. 
Most importantly, however, the systematic performance 
of our approximation(s) in Fig. 1 demonstrates that the 
magnetic-field effects, electron currents, and spin are cor- 
rectly accounted for. 

We note that the good accuracy of the LSDA in terms 
of total energies (see the dotted line in Fig. 1) is due to 
the compensation of respective errors in the exchange and 
correlation energies. 25 On the other hand, we tested our 
correlation functional for the important limit of the ho- 
mogeneous 2D electron gas (2DEG), for which the LSDA 
correlation is exact, and found reasonable agreement as 
a function of r s for both zero and full spin polarization 
(see Fig. 2). Here we used the original average parameter 
values c aa = 1.32 and c aS — 0.75. Note that according to 
the numerically exact results in Ref. 6, the ground-state 
of the 2DEG is unpolarized for < r s < 26. 

Finally, we point out, that in principle a given func- 
tional should be evaluated with KS orbitals obtained 
from self-consistent calculation instead of a post-hoc 
manner as we have done in this work. However, the vari- 
ational nature of DFT implies that if one evaluates the 
total energy with densities which slightly differ from the 
self-consistent one, the resulting change in the energy is 
of second order in the deviation of the densities. 




5 10 15 20 

r s (a.u.) 

FIG. 2: (color online). Correlation energy per electron in a 
homogeneous two-dimensional electron gas of full spin polar- 
ization (upper curves) and zero polarization (lower curves). 
The solid lines show the result from our functional with the 
original average parameter values c acr = 1.32 and c a a = 0.75. 
The dashed lines show to the local spin-density approximation 
for the correlation, corresponding to the numerically exact re- 
sult in this system. 6 



IV. CONCLUSIONS 

We have derived a spin- and current-dependent ap- 
proximation for the correlation energy of finite two- 
dimensional electron systems. The core of the deriva- 
tion is a model for the correlation-hole function of both 
same-spin and opposite spin pairs, respectively, that sat- 
isfies a set of exact properties. The excellent results ob- 
tained for few-electron quantum dots with different spin- 
polarization, current, external magnetic field, and cov- 
ering a wide range of correlation energies, strongly rec- 
ommend further developments along the construction we 
have presented here. 
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